# -*- coding: utf-8 -*-
# @Author  : longbhu
# @Time    : 2025/3/24 14:45
# @Function:
import rasterio
import numpy as np

def read_tif(file_path):
    with rasterio.open(file_path) as src:
        data = src.read(1).astype(np.float32)
        meta = src.meta
        nodata_value = src.nodata
        if nodata_value is not None:
            data[data == nodata_value] = np.nan
        return data, meta

def calculate_SCF(Cl, OM):
    # 确保所有输入都是有效的（非 nan）
    valid_mask = ~np.isnan(Cl) & ~np.isnan(OM)

    SCF = np.full_like(Cl, np.nan, dtype=np.float32)
    SCF[valid_mask] = 1 / (1 + 0.0066 * Cl[valid_mask]**2 + 0.021 * OM[valid_mask]**2)

    return SCF

def save_as_tif(data, meta, output_path):
    meta.update(
        driver='GTiff',
        count=1,
        dtype=data.dtype,
        nodata=-9999
    )

    with rasterio.open(output_path, 'w', **meta) as dst:
        dst.write(data, 1)

# 读取tif文件中的变量数据及其元数据
path_to_Cl = r'F:\code\dev\gep-calculation-helper\SF\input\SR_SF_inpu\T_clay_china.tif'
path_to_C = r'F:\code\dev\gep-calculation-helper\SF\input\SR_SF_inpu\ORG_CARBON1_china_1km.tif'

Cl, meta_Cl = read_tif(path_to_Cl)
C, meta_C = read_tif(path_to_C)
OM = 1.724 * C

# 确保所有变量具有相同的形状和元数据
assert Cl.shape == OM.shape, "All input rasters must have the same dimensions"
meta = meta_Cl  # 使用其中一个变量的元数据作为输出文件的元数据

# 计算SCF值
SCF = calculate_SCF(Cl, OM)

# 将结果保存为tif文件
output_path = './output/SF_SCF_china_1km.tif'
save_as_tif(SCF, meta, output_path)

print(f"SCF value saved to {output_path}")
